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ing  completed  the  electronic  part  of  the  calculation,  one  then  solves  the  nuclear 
part  of  the  problem  either  classically  or  quantum  mechanically.  As  an  example 
of  the  sequential  ab  initio  approach,  the  infrared  and  Raman  rotational  and 
vibrational-rotational  spectral  band  contours  for  the  water  molecule  are  com¬ 
puted  in  the  simplest  rigid  rotor,  normal  mode  approximation.  Quantum  tech¬ 
niques,  are  used  to  calculate  the  necessary  potential  energy,  dipole  moment, 
and  polarizability  information  at  the  equilibrium  geometry.  A  new  quick,  accu¬ 
rate,  and  easy  to  program  classical  technique  involving  no  reference  to  Euler 
angles  or  special  functions  is  developed  to  compute  the  infrared  and  Raman 
band  contours  for  any  rigid  rotor,  including  asymmetric  tops.  A  second,  or 
simultaneous,  type  of  ab  initio  approach  is  suggested  for  large  systems,  particu¬ 
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in  which  the  electronic  and  nuclear  parts  of  the  Bom-Oppenheimer  approxima¬ 
tion  are  simultaneously  solved  is  needed.  A  quantum  force  classical  trajectory 
(QFCT)  molecular  dynamic  method,  based  on  linear  response  theory,  is 
described,  in  which  the  forces,  dipole  moment,  and  polarizability  are  computed 
quantum  mechanically,  using  gradient  techniques  step  by  step  along  a  classical 
trajectory  whose  path  is  determined  by  these  quantum  forces.  We  believe  the 
QPCT  method  to  be  a  more  practical  ab  initio  route  to  spectral  band  contours 
for  large  molecules,  dusters,  and  solutions,  and  it  can  be  equally  applied  to 
equilibrium  and  non-equilibrium  systems.  It  is  pointed  out  that  a  similar  ab  ini¬ 
tio  QFCT  molecular  dynamic  approach  could  be  used  to  compute  other  types  of 
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transport  properties  and  thermodynamic  function.'  ^nd  their  quantum  correc¬ 
tions.  For  parameters  not  depending  on  momenta,  a  parallel  ab  initio  Monte 
Carlo  approach  would  use  electronic  energies  and  other  parameters  of  interest 
generated  quantum  mechanically,  and  "classical"  trial  moves  of  the  nuclei. 
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ABSTRACT 

We  discuss  several  ways  in  which  molecular  absorption  and  scattering 
spectra  can  be  computed  ab  initio,  from  the  fundamental  constants  of  nature. 
These  methods  can  be  divided  into  two  general  categories.  In  the  first,  or 
sequential,  type  of  approach,  one  first  solves  the  electronic  part  of  the 
Schrddinger  equation  in  the  Bom-Oppenheimer  approximation,  mapping  out 
the  potential  energy,  dipole  moment  vector  (for  infrared  absorption)  and  polari¬ 
zability  tensor  (for  Raman  scattering)  as  functions  of  nuclear  coordinates.  Hav¬ 
ing  completed  the  electronic  pan  of  the  calculation,  one  then  solves  the  nuclear 
part  of  the  problem  either  classically  or  quantum  mechanically.  As  an  example 
of  the  sequential  ab  initio  approach,  the  infrared  and  Raman  rotational  and 
vibrational-rotational  spectral  band  contours  for  the  water  molecule  are  com¬ 
puted  in  the  simplest  rigid  rotor,  normal  mode  approximation.  Quantum  tech¬ 
niques,  are  used  to  calculate  the  necessary  potential  energy,  dipole  moment, 
and  polarizability  information  at  the  equilibrium  geometry.  A  new  quick,  accu¬ 
rate,  and  easy  to  program  classical  technique  involving  no  reference  to  Euler 
angles  or  special  functions  is  developed  to  compute  the  infrared  and  Raman 
band  contours  for  any  rigid  rotor,  including  asymmetric  tops.  A  second,  or 
simultaneous,  type  of  ab  initio  approach  is  suggested  for  large  systems,  particu¬ 
larly  those  for  which  normal  mode  analysis  is  inappropriate,  such  as  liquids, 
clusters,  or  floppy  molecules.  Then  the  curse  of  dimensionality  prevents  map¬ 
ping  out  in  advance  the  complete  potential,  dipole  moment,  and  polarizability 
functions  over  the  whole  space  of  nuclear  positions  of  all  atoms,  and  a  solution 
in  which  the  electronic  and  nuclear  parts  of  the  Bom-Oppenheimer  approxima¬ 
tion  are  simultaneously  solved  is  needed.  A  quantum  force  classical  trajectory 
(QFCT)  molecular  dynamic  method,  based  on  linear  response  theory,  is 
described,  in  which  the  forces,  dipole  moment,  and  polarizability  are  computed 
quantum  mechanically,  using  gradient  techniques  step  by  step  along  a  classical 
trajectory  whose  path  is  determined  by  these  quantum  forces.  We  believe  the 
QFCT  method  to  be  a  more  practical  ab  initio  route  to  spectral  band  contours 
for  large  molecules,  clusters,  and  solutions,  and  it  can  be  equally  applied  to 


equilibrium  and  non-equilibrium  systems.  It  is  pointed  out  that  a  similar  ab  ini¬ 
tio  QFCT  molecular  dynamic  approach  could  be  used  to  compute  other  types  of 
spectra,  for  example  electronic  absorption,  as  well  as  other  parameters  such  as 
transport  properties  and  thermodynamic  functions  and  their  quantum  correc¬ 
tions.  For  parameters  not  depending  on  momenta,  a  parallel  ab  initio  Monte 
Carlo  approach  would  use  electronic  energies  and  other  parameters  of  interest 
generated  quantum  mechanically,  and  "classical”  trial  moves  of  the  nuclei. 
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I.  INTRODUCTION 

With  increased  computer  power  and  improved  computational  techniques,  such  as  the  gra¬ 
dient  methods  described  in  Sections  III  and  IV  below,  it  is  becoming  practical  to  compute  spec¬ 
tra  ab  initio,  from  the  fundamental  constants  of  nature,  for  systems  of  increasing  complexity. 
Such  an  approach  is  appealing  not  only  because  of  the  usual  academic  desire  to  elegantly  derive 
observed  phenomena  from  first  principles,  but  also  because  of  the  practical  desire  to  be  able  to 
experimentally  identify  and  understand  transient  species  or  states  (such  as  those  existing  during 
the  course  of  chemical  reactions)  whose  spectral  properties  cannot  be  deduced  from  equilibrium 
measurements.  Therefore  in  this  paper  we  explore  several  possible  ab  initio  approaches  to  spec¬ 
tra,  specifically  focusing  on  infrared  and  nonresonance  Raman.  Many  of  the  points  discussed 
can,  however,  also  be  extended  to  other  types  of  spectra  such  as  electronic  absorption  and  reso¬ 
nance  Raman. 

We  discuss  two  approaches,  sequential  and  simultaneous.  The  sequential  approach,  in  which 
first  the  electronic  part  and  then  later  the  nuclear  part  of  the  Bom-Oppenheimer  approximation 
is  solved,  is  appropriate  for  small  systems.  The  simultaneous  approach,  in  which  the  electronic 
and  nuclear  parts  are  solved  at  the  same  time,  is  more  appropriate  for  many  atom  systems. 
Then  we  review  the  newer  quantum  gradient  techniques  which  greatly  simplify  ab  initio  calcula¬ 
tions  of  spectra.  As  a  simple  illustration  of  a  sequential  ab  initio  calculation,  we  compute  the 
infrared  and  Raman  spectral  band  contours  for  the  water  molecule  using  a  rigid  rotor,  normal 
mode  approach.  The  vibrations  are  handled  quantum  mechanically,  but  in  order  to  make  the 
connection  to  liquid  state  spectra  in  which  the  individual  ro vibrational  lines  have  usually 
merged,  we  compute  here  the  band  contours.  While  this  could  be  done,  as  illustrated  else¬ 
where,1*2  by  broadening  the  quantum  rovibrational  lines,3*4  we  choose  instead  in  this  paper  to 
explore  another  path,  computing  the  band  contours  classically  in  the  rigid  rotor  limit,  with  the 
aim  of  eventually  gaining  a  more  intuitive  understanding  of  the  relation  among  moments  of 
inertia,  dipole  and  polarizability  derivatives,  temperature,  and  band  contour  shape.  The  classical 
technique  we  have  developed  is  a  simple,  quick,  and  easy  way  to  compute  band  contours,  and 
might  be  used,  for  example,  to  conveniently  explore  the  dipole  and  polarizability  components 
consistent  with  measured  spectra.  Once  the  rotational  correlation  functions  have  been 
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computed  from  the  moments  of  inertia  and  Fourier  transformed,  one  could  compute  band  con¬ 
tours  immediately  for  any  combination  of  dipole  and  polarizability  derivatives  simply  by  weight¬ 
ing  these  transformed  correlation  functions  by  the  appropriate  dipole  or  polarizability  com¬ 
ponents  and  summing.  An  on-line  visual  matching  of  theoretical  to  experimental  infrared  and 
Raman  band  contours  is  quite  possible  and  might  be  very  instructive  in  understanding  the  phy¬ 
sical  basis  of  observed  spectra,  for  example  the  symmetries  and  directions  of  dipole  moment 
and  polarizability  changes  associated  with  different  vibrations.  This  can  be  carried  out  even  in 
the  absence  of  potential  energy  surface  information  of  sufficient  accuracy  to  match  line  spectra. 
In  this  way,  for  example,  one  could  use  crude  ab  initio  dipole  and  polarizability  parameters  to 
help  resolve  the  usual  experimental  sign  uncertainties  and  then  refine  the  parameters  by  more 
detailed  comparison  with  experimental  band  contours.  We  derive  the  appropriate  equations  to 
compute  infrared  and  Raman  spectra  and  present  our  ab  initio  sequential  rigid  rotor,  normal 
mode  bond  contour  results  for  Hfi  and  compare  with  experiment.  Finally  we  describe  in  more 
detail  our  proposed  ab  initio  simultaneous  approach,  which  can  in  principle  allow  the  computation 
of  spectra  and  many  other  properties  for  larger  systems  using  a  dose  linking  of  quantum 
evaluation  of  forces  and  other  parameters  with  classical  trajectory  integration. 

II.  SEQUENTIAL  VERSUS  SIMULTANEOUS  APPROACHES 

In  the  usual  sequential  approach,  the  electronic  part  of  the  Bom-Oppenheimer  approxima¬ 
tion  is  first  solved  for  the  potential  energy  V  as  a  function  of  the  positions  of  the  nuclei,  and 
for  the  connection  to  the  radiation  field,  for  example  the  dipole  moment  vector  /*.  for  infrared 
spectra  and  the  polarizability  tensor  P  for  Raman  spectra.  The  least  complicated  sequential  ab 
initio  approach,  which  we  illustrate  here  by  a  calculation  of  the  infrared  and  Raman  spectral 
band  contours  for  the  water  molecule,  is  to  search  for  the  minimum  of  the  potential  energy 
surface  and  then  to  expand  V,  p,  and  P  in  Taylor’s  series  about  this  minimum.  At  the  sim¬ 
plest  level,  a  rigid  rotor,  normal  mode  approximation  can  be  made  for  the  nudear  rotational  • 
vibrational  (rovibrational)  wavefunctions  and  the  derivatives  of  ja  and  P  with  respect  to  the 
normal  coordinates  can  then  be  used  to  compute  the  intensities  of  the  various  rovibrational 
transitions.1'2’ 5-6  Higher  derivatives  in  the  Tabor’s  series  for  V,  jt,  and  P  can  be  used  in  more 
accurate  calculations.  An  even  more  accurate  sequential  solution,  applicable  if  the  molecule  is 
quite  small,  is  to  solve  for  V,  jt,  and  P  over  all  the  necessary  coordinate  space  and  then 
numerically  to  solve  for  the  rovibrational  eigenfunctions  and  energy  eigenvalues,  and  finally 
numerically  to  solve  for  the  transition  intensities.  An  interesting  and  instructive  discussion  of 
various  theoretical  approaches  to  the  computation  of  infrared  spectra  with  particular  application 
to  water  monomers,  dimers,  and  liquid  water  has  recently  been  given  by  Coker,  Reimers,  and 
Watts.7  An  intriguing  ab  initio  calculation  for  the  N2  molecule  has  been  carried  out  by  Svendsen 
and  Oddershede*  who  use  the  frequency  dependent  polarizability  tensor  instead  of  the  usual 
static  one.  Several  other  types  of  semidassical  sequential  approaches  could  also  be  carried  out  in 
an  ab  initio  fashion.9' 10 

For  molecular  systems  of  more  than  a  few  atoms,  the  simultaneous  approach  to  ab  initio 
calculations  is  preferred.  As  the  size  of  the  system  increases,  it  becomes  increasingly  impracti¬ 
cal  to  link  quantum  mechanics  to  experimental  observables  through  the  complete  multidimen¬ 
sional  computation  and  storage  of  intermediate  functions  of  nuclear  coordinates  such  as  poten¬ 
tial  energy  V,  dipole  moment  p,  and  polarizability  P.  As  the  number  N  of  atoms  grows,  the 
"curse  of  dimensionality”  makes  the  explicit  computation  as  well  as  the  fitting  and  storage  of 
such  functions  infeasible.  For  example,  if  10  grid  points  are  needed  for  sufficient  accuracy  in 
each  dimension,  then  l(fiN  evaluations  are  needed  to  specify  the  full  function,  or  10*  if  a 
dozen  atoms  are  involved.  Clearly,  no  experimental  measurement  contains  or  can  require  this 
much  information,  and  a  different  theoretical  approach  is  called  for.  An  alternative  is  statistical 
sampling,  for  example  using  a  molecular  dynamics  or  Monte  Carlo  approach.  Specifically,  as 
one  of  us  has  suggested,11  classical  molecular  dynamics  may  be  integrated  with  ab  initio  quan¬ 
tum  force  evaluation,  so  that  at  each  time  step  the  force  on  each  nucleus  is  evaluated  quantum 
mechanically  and  then  used  to  integrate  forward  one  step  at  a  time  the  classical  trajectory.  One 

would  follow  the  trajectory  in  time,  sampling  the  force  F,  -  ~jr— ,  the  dipole  moment  p,  and 


the  polarizability  P  only  at  specific  nuclear  configurations  through  which  the  trajectory  passes, 
and  then  apply  the  techniques  illustrated  in  previous  papers1*2  to  derive  the  infrared  and 
Raman  spectra  from  the  power  spectra  of  the  time  histories  of  jt(r)  and  P(r).  While  this  sam¬ 
pling  may  build  up  over  the  course  of  many  trajectories  a  continuously  improving  estimate  of 
the  entire  V,  n  and  P  functions  over  the  accessible  volume  of  configuration  space,  we  do  not 
initially  require  a  global  representation  of  the  functions,  and  in  systems  of  many  atoms  we  are 
likely  to  achieve  sufficient  accuracy  to  compare  with  experiment  long  before  such  a  global 
representation  could  be  computed.  With  an  average  over  the  proper  ensemble  of  initial  nuclear 
positions  and  momenta,  observables  ranging  from  thermodynamic  functions  to  spectra  to  chem¬ 
ical  reaction  properties  to  transport  properties  could  be  computed  in  this  way.  An  alternative 
approach,  for  phenomena  which  are  not  explicitly  time  dependent,  is  a  Monte  Carlo  procedure 
in  which  at  each  trial  nuclear  configuration  step  the  energy  is  evaluated  quantum  mechanically, 
and  for  accepted  steps  the  other  variables  of  interest  are  then  computed. 

Two  initial  calculations  of  this  type,  which  for  the  molecular  dynamics  case  we  will  call 
Quantum  Force  Classical  Trajectory  (QFCT),  have  been  carried  out  by  Warshd  and  Karplus12 
and  by  Leforestier.13  Warshel  and  Karplus12  use  forces  given  from  analytical  first  derivatives  in 
a  semidassical  technique  to  integrate  the  equations  of  motion  for  cis-tnms  photoisomerization 
from  the  triplet  state  of  butene-2.  Leforestier13  studied  the  gas  phase  nucleophilic  substitution 
reaction  H~  +  CH4  —CH4  +  H~.  He  derived  his  forces  using  an  ab  initio  approach  which  will 
be  reviewed  below  and  integrated  the  classical  trajectories  outward  from  an  initial  state 
corresponding  to  the  transition  state  for  the  reaction. 

Thus,  classical  statistical  mechanics  plus  an  ab  initio  QFCT  method  can  be  used  to  predict 
spectral  observables,  computing  the  time  evolution  of  the  dipole  moment  vector  jt(r)  and  the 
polarizability  tensor  P(/),  and,  as  we  have  previously  demonstrated,  from  ensemble  averages  of 
their  power  spectra  the  infrared1  and  Raman2  spectra.  This  can  be  extended  to  computing 
infrared  and  Raman  spectra  for  more  complex  systems,  such  as  clusters  and  liquid  solutions, 
and  transient  spectra  during  chemical  reactions,14  or  in  other  words  to  situations  for  which  only 
bond  contours  may  be  observable  experimentally  and  for  which  trying  to  quantize  nuclear 
motion  would  be  too  difficult  theoretically.  We  emphasize  i)  that  other  spectral  and  transport 
phenomena  could  also  similarly  be  derived  through  calculation  of  their  proper  correlation 
functions,15'1*  ii)  that  thermodynamic  variables  and  their  quantum  corrections  could  be  com¬ 
puted,19  and  iii)  that  the  parallel  Monte  Carlo  approach  for  properties  independent  of  momenta 
would  be  quantum  evaluation  of  energy  (plus  force20'22  if  desired)  at  each  trial  step  followed  by 
evaluation  of  other  variables  to  be  averaged  for  those  configurations  which  are  accepted. 

III.  QUANTUM  POTENTIAL  ENERGY,  FORCE,  AND  EQUILIBRIUM  GEOMETRY 

The  emergence  of  gradient  methods  in  recent  years  allows  the  use  of  very  efficient 
methods  for  the  calculation  of  forces,  the  location  of  equilibrium  structures  and  the  evaluation 
of  harmonic  force  constants  and  vibrational  frequencies.23*30  More  recently,  one  of  us  has  sug¬ 
gested  that  the  power  of  these  methods  also  be  used  to  evaluate  the  dipole  and  polarizability 
derivatives.31'32  (For  an  example  of  infrared  vibrational  properties  computed  without  gradient 
techniques  see  Person,  Brown,  Steele,  and  Peters.33*34  )  In  Sections  III  and  IV  we  summarize 
the  salient  points  of  the  gradient  approach  and  show  how  these  methods  are  used  to  evaluate 
the  relevant  quantities  required  in  the  example  sequential  ab  initio  spectral  analysis  for  the  H2O 
infrared  and  Raman  band  contours  presented  below.  The  molecular  calculations  reported  here 
are  all  performed  at  the  SCF  level  with  the  GRADSCF  program  system.35 

In  many  chemical  applications  one  needs  to  evaluate  a  molecular  wavef unction  from 
which  the  energy,  dipole  moment  and  other  relevant  properties  are  calculated.  It  was  recognized 
a  number  of  years  ago36*39  that  often  a  knowledge  not  only  of  the  energy  U,  but  also  of  its  first 
derivatives  with  respect  to  the  nuclear  coordinates  can  considerably  enhance  the  scope  of  chem¬ 
ical  computational  investigations.  Over  the  past  decade  efficient  means  of  evaluating  the  first 
derivatives  of  the  energy  have  been  developed,  first  for  semiempirical24  and  eventually  for  ab 
initio  wavef unctions. 231 28-  40  In  these  applications  both  the  energy  and  the  gradient  are  evaluated 
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simultaneously.  The  derivation  of  the  gradient  expression  Tor  a  variety  of  wavefunctions  has 
been  extensively  discussed  in  the  literature23- 2I- 40  and  will  not  be  repeated  here. 

If  we  expand  the  energy  as  a  function  of  the  nuclear  coordinates  in  a  Taylor's  series  we 
find  that 

UCX)  -  U(Xo)  +  gf(X  -  Xo)  +  (X  -  XoTQCX  -  Xo)  (3.1) 

where  X  is  the  vector  of  Cartesian  nuclear  coordinates,  g  is  the  gradient  vector  whose  elements 

Aft 

are  defined  as  gj  —  -™,  and  Q  is  the  Hessian,  or  matrix  of  second  derivatives.  The  geometry 

oXj 

optimization  is  performed  by  employing  the  following  algorithm, 

X-Xo-oQ-**,  (32) 

where  a  is  a  scale  factor.  The  specific  minimization  algorithm  which  we  use  is  that  of  Mur- 
taugh  and  Sargent  41  where  the  matrix  Q  is  initially  approximated  by  a  unit  matrix  and  updated 
after  each  successive  step  to  form  the  inverse  to  the  true  Hessian.  Convergence  in  the 
geometry  optimization  is  deemed  satisfactory  if  the  largest  component  of  the  gradient  is  less 
than  4x1 0-4  hartree/bohr. 

The  force  constants  are  evaluated  by  the  finite  difference  techniques  described 
previously.25'27  By  numerically  differencing  the  gradient  vector  we  are  able  to  evaluate  a 
column  of  the  matrix  of  second  derivatives  for  each  displacement  of  the  nuclear  coordinates. 
The  force  constants  are  evaluated  in  a  Cartesian  representation  and  the  frequencies  and  normal 
modes  are  determined  by  diagonalizing  the  mass  weighted  force  constant  matrix.5-25*27  We  also 
evaluate  the  appropriate  B  matrix5-6-25  which  in  turn  is  used  to  transform  the  force  constants 
from  a  Cartesian  to  an  internal  symmetry  coordinate  representation.  This  transformation  serves 
rigorously  to  project  out  any  translation  and  rotation  from  the  vibrational  modes  and  yields  six 
true  zero  frequency  modes. 

The  water  molecule  results  for  the  equilibrium  geometry,  moments  of  inertia,  and  vibra¬ 
tional  frequencies  derived  from  harmonic  force  constants  are  shown  in  Table  I.  Here  r  is  the 
O-H  distance,  9  is  the  H-O-H  angle,  Iatt  is  the  aa  component  of  the  moment  of  inertia  tensor, 
A  is  Angstrom,  deg  is  degrees,  amu  is  atomic  mass  units,  and  cm-1  is  wavenumbers  in  inverse 
centimeters.  Figure  1  shows  the  Cartesian  axis  system  chosen  for  the  water  molecule,42  and 
Fig.  2  shows  the  normal  modes.42  The  symmetric  stretch  corresponds  to  v\  and  qi,  v2  and  q2 
correspond  to  the  symmetric  bend,  and  and  qj  correspond  to  the  asymmetric  stretch,  in 
which  Vi  denotes  the  ith  normal  mode  frequency  and  qj  is  the  ith  normal  coordinate. 


TABLE  I.  Ab  Initio  parameters  computed  for  HjO  molecule. 


r 

9 

4r 

4 

4 

v\ 

*2 

vj 

(A) 

(det) 

( amu  A2) 

(amu  A2) 

(i amu  A2) 

(cm-1) 

(cm-1) 

(cm-1) 

0.943 

106.0 

1.720 

0.5766 

1. 1434 

3797 

1740 

3902 

IV.  QUANTUM  DIPOLE  MOMENT  AND  POLARIZABILITY  TENSOR 

The  dipole  and  polarizability  derivatives  with  respect  to  the  normal  modes  are  computed 
using  a  very  efficient  new  method  which  has  been  reported  previously.31-32  We  seek  for  the 
infrared  spectrum  the  quantities 


where  pa,  a  -  x,y,z,  are  the  three  components  of  the  molecule  fixed  dipole  moment  vector  j* 
and  the  are  the  normal  mode  coordinates.  For  the  polarizability  derivatives  we  seek  the 
analogous  quantities 
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(4.2) 


where  the  Pap  are  the  components  of  the  molecule  fixed  polarizability  tensor,  in  which  a  and 
/9  —  x,y,z  and  again  the  q  are  the  normal  coordinates.  The  dipole  moment  vector  and  the 
polarizability  tensor  are  the  first  and  second  derivatives  of  the  energy  U(E)  in  the  presence  of 
an  applied  external  electric  field  E.  We  can  therefore  write  the  derivatives  of  the  elements  of 
the  dipole  moment  vector  and  polarizability  tensor  with  respect  to  the  Cartesian  components  Xj 

ai  t 

of  the  positions  of  the  nuclei  in  the  following  manner,  in  which  as  before  &  —  , 


3m»  8  8U(E) 8  8U(E) 8_  ,-v 

9Xj  8Xi  8E.  8Ea  8Xj  “  8E„  * 


3P.0  8  aHJ(E)  a2  8U(E)  82  fn 

8Xj  “  9Xj  8Ea8Efl  “  8Ea9E„  9Xj  "  3Ea8E/l 


(4.3) 

(4.4) 


The  final  terms  in  both  equations  result  from  the  interchange  in  the  order  of 
differentiation.  These  final  expressions  involve  the  derivatives  of  the  gradient  of  the  energy 
with  respect  to  the  components  of  an  applied  electric  field.  If  we  now  transform  the  resulting 
gradient  vector  from  a  Cartesian  to  a  normal  coordinate  representation  we  arrive  at  the  desired 
quantities  shown  in  Eqs.  (4.1)  and  (4.2).  The  simplicity  and  efficiency  of  the  method  lies  in 
the  fact  that  the  introduction  of  the  electric  field  involves  only  a  one  electron  perturbation  to 
the  molecular  Hamiltonian.  We  write  the  Hamiltonian  as 

H(E)  -  Ho  +  H'(E)  (4.5) 

where 

H'(E)  -  e^ETj  -  eY Zj  E  Rj  (4.6) 

i-i  j-i 

in  which  e  is  the  magnitude  of  the  electron  charge,  E  is  the  applied  electric  field,  r{  is  the  posi¬ 
tion  of  the  ith  electron,  the  first  sum  is  over  the  electrons,  Zj  is  the  charge  of  the  jth  nucleus, 
Rj  is  its  position,  and  the  second  sum  is  over  the  nuclei.  For  the  computation  of  the  energy 
this  Hamiltonian  requires  the  addition  of  a  term  which  includes  the  appropriate  dipole  moment 
integrals.  It  is  important  to  realize  that  the  variation  of  the  field  is  independent  of  the  nuclear 
coordinates.  In  all  of  our  calculations  for  the  water  molecule  presented  here  the  molecule  is 
held  fixed  at  the  equilibrium  geometry  and  SCF  calculations  are  performed  for  a  set  of  applied 
electric  fields.  For  any  polyatomic  system,  the  dipole  derivatives  require  at  most  six  SCF  calcu¬ 
lations:  two  variations  of  the  applied  field  in  each  of  the  three  Cartesian  directions.  The  polari¬ 
zability  derivatives  require  a  total  of  12  SCF  calculations,  since  these  require  a  two  dimensional 
difference  formula  for  the  finite  difference  equations.  All  of  these  calculations  use  the  same  set 
of  one  and  two  electron  integrals.  Furthermore,  the  evaluation  of  the  gradient  is  performed 
only  once.  In  the  absence  of  an  applied  field  the  gradient  calculation  yields  one  vector.  Once  we 
apply  the  field  in  all  of  the  required  directions,  the  gradient  calculation  yields  either  six  or 
twelve  vectors. 

These  gradient  vectors,  which  are  now  a  function  of  the  electric  field,  are  then 
transformed  from  a  Cartesian  to  a  normal  coordinate  representation.  The  resulting  values  are 
numerically  differentiated  to  yield  the  dipole  and  polarizability  derivatives. 

In  the  geometry  and  force  constant  calculations  for  the  water  molecule  we  use  a  split- 
valence  polarized  basis  which  was  originally  developed  by  Pople  and  co-workers.43  This  basis  set 
is  described  as  a  (10s4pld/3slp)  contracted  to  a  [3s2pld/2slp]  basis  and  often  referred  to  as  a 
6-3IG**  basis.  The  geometries  and  vibrational  frequencies  calculated  with  this  basis  are  close  to 
the  limit  attainable  at  the  SCF  level.  The  d-type  polarization  functions  on  the  oxygen  atom  have 
an  exponent  of  0.8,  while  the  p-type  polarization  functions  on  the  hydrogen  atoms  are  given 


exponents  of  1.1. 

The  calculation  of  the  water  dipole  and  polarizability  derivatives  requires  a  much  more 
extensive  basis  set.  It  has  previously  been  shown  that  a  targe  and  flexible  basis  set  is  needed  to 
describe  accurately  the  interaction  of  a  molecule  with  an  applied  electric  fleld.  To  this  end  we 
employ  a  triple-zeta  plus  double  polarization  basis  (TZPP)  for  the  dipole  and  polarizability 
derivatives.  This  basis  set  provides  a  near  Hartree-Fock  description  of  the  molecular  bonding, 
plus  a  set  of  diffuse  functions  of  sufficient  flexibility  to  interact  with  the  applied  field.  In  these 
calculations  we  employ  a  Huzinaga  basis  of  (lls6p/5s)  as  contracted  by  Dunning44  to  a 
I5s3p/3s]  basis  on  the  oxygen  and  the  hydrogen  atoms.  This  basis  is  augmented  by  two  sets  of 
d  functions  on  the  oxygen  with  exponents  of  0.8  and  0.078,  and  two  sets  of  p-type  polarization 
functions  on  the  hydrogens  with  exponents  of  1.1  and  0.116.  The  magnitude  of  the  applied 
field  is  chosen  to  be  l.OxlO-3  au.  The  results  of  the  calculation  are  presented  in  Tables  II  and 
III,  in  which  na  is  the  a  component  of  the  dipole  moment  vector,  Paa  is  the  aa  diagonal  ele¬ 
ment  of  the  polarizability  tensor  P,  D  is  debye,  and  A  is  angstrom.  Extensive  comparison  with 
experimentally  derived  values  of  these  parameters  will  not  be  presented  in  this  paper,  and  we 
only  briefly  note  here  that  the  dipole  parameters  may  be  compared  with  the  experimental  values 
of  Camy-Peyret  and  Flaud45  and  the  polarizability  parameters  with  the  work  of  Murphy.46-4* 

TABLE  II.  Non-zero  ab  initio  equilibrium  dipole  moment 

and  polarizability  components. 


TABLE  III.  Ab  initio  dipole  and  polarizability  derivatives. 
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2> 

(D  A"’  amu 
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1.4663 

0.5538 

t.5355 

0.0 

1 

(k2amu  2) 
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-I 

(k2amu  2) 

0.40349  0.0 

0.0 

0.0252 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

2.03936 

0.0 

0.0 

0.18969 

0.0 

0.0 

0.0 

1.1048 

0.0 

0.0 

1.3186 

0.0 

0.0 

-0.29568 

0.0 

1.1048 

0.0 

V.  CLASSICAL  ROTATIONAL  TRAJECTORIES 

Given  the  moments  of  inertia  computed  from  the  equilibrium  coordinates  of  the  water 
molecule  as  shown  in  Section  II!  above,  we  could  straightforwardly  compute  the  rotational  spec¬ 
trum3- 6  quantum  mechanically  in  the  rigid  rotor  approximation.  Instead,  in  order  to  gain 
more  intuitive  insight  into  the  physical  basis  of  infrared  and  Raman  band  contours,  and  to 
make  more  evident  the  connection  to  the  liquid  state,  we  compute  the  band  contours  from  clas¬ 
sical  mechanics,  using  a  method  improved  over  the  asymmetric  rotor  techniques  pioneered  by 
Bratos,  Guissani,  and  Leicknam. 49-56  The  methodology  presented  in  Sections  V-VII  allows  the 
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rapid  and  accurate  calculation  of  classical  rigid  rotor  infrared  and  Raman  band  contours  for  any 
molecule,  including  asymmetric  rotors  like  water.  (For  examples  of  earlier  less  general  treat¬ 
ments  for  linear,  spherical,  and  symmetric  tops,  see  St.  Pierre  and  Steele57  and  references 
therein  as  well  as  references  in  the  asymmetric  rotor  papers  of  Bratos,  Guissani,  and 
Leicknam.49*56  )  The  method  presented  here  is  simple  and  direct,  involving  no  reference  to 
Euler  angles  or  special  functions.  Approximate  sets  of  infrared  and  Raman  contours  can  be 
computed  in  one  minute  on  a  VAX58  11/780  and  very  accurate  sets  in  IS  minutes.  The  rota¬ 
tional  functions  discussed  here  require  only  a  250  line  program  plus  standard  integration  and 
fast  Fourier  transform  routines.  In  the  rest  of  this  section  we  develop  the  techniques  to 
efficiently  compute  the  necessary  rotational  properties. 

In  order  to  calculate  rotational  trajectories  we  choose  a  set  of  body  fixed  axes  coinciding 
with  the  principle  axes  of  the  molecule,  those  axes  of  a  rigid  body  for  which  the  moment  of 
inertia  tensor  is  diagonal.  The  center  of  mass  is  at  the  origin  of  both  the  space  fixed  and  body 
fixed  coordinate  systems;  thus,  the  relationship  between  the  two  sets  of  axes  is  specified  by 
their  relative  orientation.  This  orientation  is  usually  expressed  in  terms  of  the  Euler  angles  <M, 
and  but  for  our  purposes  (particularly  the  numerical  aspects)  it  is  more  convenient  to  use 
the  Cartesian  rotation  matrix  D  in  terms  of  direction  cosines  (for  an  alternative  treatment  using 
Euler  angles  and  special  functions,  see  the  work  of  Bratos,  Guissani,  and  Leickman49*56  ).  If  a 
is  a  vector  in  terms  of  the  space  fixed  axes,  and  a'  is  the  vector  in  terms  of  the  body  fixed  axes, 
then  a  is  related  to  a*  by  a  -  Da'.  Similarly,  if  T  is  a  tensor,  we  have  T  -  DTD-1.  The  time 
history  of  the  orientation  of  the  molecule  is  given  by  D(r). 

We  now  derive  equations  of  motion  for  D(t).  While  the  results  are  not  new,  the  follow¬ 
ing  derivation  is  ample  and  serves  to  introduce  the  notation  used  in  later  parts  of  the  paper. 
Let  e,,  /— 1,2,3,  be  the  unit  vectors  in  the  x,  y,  and  z  space  fixed  directions,  and  similarly  let 
E„(f),  a— 1,2,3,  be  the  time  varying  body  fixed  unit  vectors.  Note  that  we  use  lower  case 
Roman  letters  for  indices  taking  on  values  associated  with  the  space  frame,  and  lower  case 
Greek  letters  for  body  frame  indices.  The  components  of  D(r),  DiaU ),  are  given  by 

A.(f) -«,•£.<*>.  (5.1) 

Let «(/)  be  the  angular  velocity  vector  and  »„(f)  be  its  components  in  the  body  frame,  so  that, 
adopting  here  and  in  the  rest  of  the  paper  (except  where  otherwise  noted),  the  convention  that 
repeated  indices  on  different  symbols  are  to  be  summed  over,  «*(f)  -  cu0(/)E0(r).  If  we  let  L 
be  the  angular  momentum,  and  I  be  the  moment  of  inertia  tensor,  then  L  -  Im,  or 

L  -  /.*»*(/)£,(/),  (5.2) 

where  Iaf  are  the  components  of  I  in  the  body  axes,  and  thus  do  not  vary  in  time  in  our  rigid 
body  approximation.  Since  for  our  isolated  system  there  are  no  external  torques, 

£  “  + - 0  (s-3> 


Since  the  t„  are  fixed  in  the  body  frame, 

— ^  “  -  UpCpayty  .  (5.4) 

Here  « ^  is  the  permutation  symbol  or  Levi-Civita  density,59  defined  to  be  zero  if  any  two  of 
its  indices  are  equal,  and  otherwise  +1  or  -1  according  to  whether  n ay  is  an  even  or  odd  per¬ 
mutation  respectively  of  1,2,3.  Using  this  relationship  and  Eq.  (5.3)  we  obtain 


0. 


(5.5) 


For  our  chosen  set  of  body  axes  the  inertia  tensor  is  diagonal,  so  this  reduces  to  Euler’s  equa¬ 
tions  of  motion  for  a  rigid  body  free  of  external  torques: 

dlttx 

I XX  (/—  lyy)  “  0, 
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-4>-0,  (5.6) 

+  •»*••>(/*  ~  4r )  “  0. 

We  similarly  obtain  from  Eqs.  (5.1)  and  (5.4) 

^A.-*ryt.-«*«#A,  (5-7) 

as  the  equations  of  motion  for  A«- 

In  order  to  calculate  infrared  and  Raman  rotational  and  vibrational-rotational  spectral  band 
contours  with  this  classical  rigid  rotor  method,  we  will  show  in  Sections  VI  and  VII  below  that 
we  must  separate  an  ensemble  average  of  an  autocorrelation  function  of  the  dipole  moment  jt 
for  the  infrared  case  and  the  polarizability  tensor  P  for  the  Raman  case  into  a  product  of 
ensemble  averages.  One  of  these  ensemble  averages  is  an  average  over  vibrational  states  of  an 
autocorrelation  function  of  or  P,  and  the  other  is  a  rotational  ensemble  average  of  elements 
of  a  rotation  matrix  over  initial  angular  velocity  and  initial  orientation.  With  Eqs.  (5.6)  and 
(5.7)  we  have  at  all  that  we  need,  in  principle,  to  calculate  numerically  the  rotational  ensemble 
averages  necessary  to  compute  the  rigid  body  rotational  part  of  the  spectra,  but  we  can  save 
ourselves  a  large  amount  of  computer  time  and  increase  our  understanding  of  the  dependence 
of  band  shape  on  the  moment  of  inertia  tensor  and  temperature  by  doing  part  of  the  problem 
analytically. 

With  Eqs.  (5.6)  and  (5.7),  D(r)  is  determined  by  the  initial  values  of  D  and  m,  i.e.  D° 
and  a»°,  while  m (/)  is  determined  solely  by  a*°.  Thus  we  write  D(r,D°,w°)  and  «*(/,m°).  If  the 
pair  Dir,!)0,*0)  and  at(/,«°)  together  satisfy  the  equations  of  motion,  then  by  the  isotropy  of 
space,  so  do  RD(r,D°,a»°)  and  «(r,»°),  where  R  is  any  rotation.  Accordingly  we  introduce  a 
rotational  Green’s  function  G(r)  with  the  following  definition: 

G(r,*°)  m  D(r,l,a*°) .  (5.8) 

in  which  1  is  the  unit  matrix.  It  follows  that  since  D(r,D°,«°)  and  D°G(r,*°)  satisfy  the  same 
equations  of  motion,  and  have  the  same  initial  conditions,  then 

D(/,D°,«*°)  -  D°G(r,*°)  .  (5.9) 

G(r)  may  be  viewed  as  a  rotational  propagator  acting  on  the  initial  orientation  D°.  G(r)  is  the 
rotation  matrix  as  a  function  of  time  assuming  that  the  initial  body  fixed  axes  coincide  with  the 
space  fixed  axes.  D°  then  rotates  this  trajectory  to  make  it  coincide  with  the  actual  initial  orien¬ 
tation.  This  separation  will  allow  us  to  analytically  perform  ensemble  averages  over  initial 
orientations,  as  we  will  show  in  Sections  VI  and  VII  below.  The  rotational  ensemble  average 
for  G  then  reduces  to  an  ensemble  average  over «°. 

As  we  will  show  below  in  Section  VI,  in  order  to  calculate  infrared  spectra  we  need  to  cal¬ 
culate  the  time  history  of  <Gm0(t)>.  Here  <  >  indicates  an  ensemble  average  over  a 
Boltzmann  distribution  of  initial  angular  velocities.  For  Raman  spectra,  as  we  will  show  in  Sec¬ 
tion  VII,  we  need  to  calculate  the  time  history  of  <  Gag(t)Grt(t)>.  Using  Eqs.  (5.6)  and 
(5.7)  we  could  calculate  these  elements  as  they  stand.  A  naive  approach  would  be  to  choose  a 
grid  in  ••  space  of  initial  w's,  calculate  the  various  elements  and  products  of  G(r)  for  a  set  of 
t’s,  weight  them  by  the  appropriate  Boltzmann  factor  and  then  sum.  (Note  that  the  equations 
of  motion  for  G  and  D  are  the  same,  since  they  differ  only  by  initial  conditions.)  Although  the 
averaging  over  initial  orientations  is  done  analytically,  the  computer  time  required  for  this  cal¬ 
culation  would  be  large,  but  fortunately  a  variety  of  simplifications  are  possible. 

Our  ensemble  average  over  initial  angular  velocity  is  given  by 

m  m  m 

<G(r)>  -  SSS  «»;.  G  (t  ,***)  dtoldutydat?  , 


(5.10) 


where  GO)  stands  for  a  term  of  the  form  Ga0(t)  for  infrared  spectra  or  Ga0(t)GrSU)  for 
Raman  spectra,  and  P(m^,m^,mf)  is  the  Boltzmann  probability  density,  which  as  we  show  in 
Appendix  A,  is  given  by 

P(m0,m°,m°)  -  Cexp(-i/3[/„K°)2  +  /„(«®)J  +  4U®)2!)  ,  (5.11) 

where  C  is  a  normalizing  constant,  /3 —{kB  T)~\  kB  is  Boltzmann’s  constant,  and  T  is  the  tem¬ 
perature.  We  can  simplify  this  probability  density  by  changing  variables  to  ua  —  -JT^uta  (no 
summation).  We  then  have 

Piu^Uy,!^)  -  C'exp|- j/J(u,J  +  u }  +  «rJ)j ,  (5.12) 

in  which  C  is  a  new  constant.  This  probability  distribution  is  now  spherically  symmetric,  so  we 
let  u1  -  u}  +  u}  +  uz2,  and  change  variables  agaiu  to  get 

P(u,Sl)  -  C"  M^xpl-I^n2! ,  (5.13) 

where  O  is  the  solid  angle,  i.e.  orientation,  of  the  vector  (ux,ufyuz).  The  ensemble  average 
then  becomes 

<G(r)>  -  Cm /wJ«"2*^/G(f,u,n)</n  du,  (5.14) 

where  GO,u,fi)  is  GO)  with  the  initial  value  m°  specified  by  u  and  O.  This  change  of  vari¬ 
ables  in  the  probability  density  means  that  we  separate  the  integration  over  the  magnitude  of  »° 
from  the  integration  over  the  direction  of  m°.  This  in  turn  will  allow  us  to  use  the  scaling  pro¬ 
perties  of  m  and  G  to  change  the  integral  over  the  magnitude  of  «  to  an  integral  over  time, 
which  will  mean  that  the  work  in  calculating  the  time  history  of  GO)  and  doing  the  integral  will 
overlap. 

We  do  this  as  follows:  we  first  note  that  solutions  to  our  equations  of  motion  are  unique 
for  given  initial  conditions.  If  we  have  two  solutions  which  have  the  same  initial  conditions, 
then  they  must  be  equal.  Using  this  property,  we  will  derive  scaling  relations  between  t  and  »°, 
first  for  at,  and  then  for  G. 

Suppose  mO,m°)  satisfies  Eqs.  (5.6).  Then  it  is  not  difficult  to  check  that  am(af,m°)  also 
satisfies  Eqs.  (5.6),  where  a  is  an  arbitrary  scalar  constant.  The  initial  conditions  for  at  tell  us 
that  aat(0,at°)  —  am0.  By  definition,  mO,am°)  is  the  solution  to  Eqs.  (5.6)  for  the  initial  value 
aat°.  We  conclude  that 

at(r,aat°)  “  am(at,m°)  .  (5.15) 

The  intuitive  picture  for  this  is  as  follows:  suppose  we  allow  two  identical  rigid  bodies  to  rotate 
subject  to  the  same  initial  conditions,  except  that  one  starts  with  an  initial  angular  velocity  in 
the  same  direction  as  the  other  but  with  a  magnitude  larger  by  a  factor  a.  This  equation  tells 
us  that  the  two  bodies  will  follow  the  same  rotational  trajectory  through  «-space,  except  that 
one  will  trace  its  path  more  quickly  by  the  factor  a  and  its  magnitude  will  be  greater  by  the  fac¬ 
tor  a.  Using  this  relation  we  can  in  a  very  similar  manner  show  that 

G(a/,m°)  -  GO, am0)  .  (5.16) 


This  relationship  tells  us  that  our  two  identical  bodies  with  different  initial  w’s  will  trace  the 
same  path  through  orientation  space  (without  scaling  the  magnitude  of  G),  except  that  one  will 
trace  its  path  more  quickly  by  the  factor  a.  We  now  change  variables  in  Eq.  (5.14)  from  u  to 


|,  with  u  -  u°(,  where  f  is  a  unitless  integration  variable,  and  where  we  choose  u0  to  be  the 
value  of  u  at  the  maximum  of  P(u,(l),  so  that  u0  —  (2//3)  .  We  obtain 


(5.17) 
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where  the  left  integral  is  over  (  from  0  to  <»  and  the  right  integral  carries  fl  over  the  sphere. 
Eq.  (S.16)  tells  us  that 

G(t,u<£,Cl)  -  <7(f/,Uo,ft)  ,  (5.18) 

so  we  have 

<G(t)>  -  *fi''ififG(et,u0,Cl)dndt ,  (5.19) 

where  we  now  explicitly  write  out  the  normalizing  constant  C*.  Written  in  this  way,  the  work 
necessary  to  compute  the  time  history  and  to  compute  the  time  integral  overlap  almost  com¬ 
pletely. 

By  various  symmetry  arguments  relating  to  the  molecule's  inertia  tensor  (which  are  valid 
because  in  the  body  frame  the  inertia  tensor  is  diagonal)  we  show  in  Appendix  B  that  the  off- 
diagonal  elements  of  <  G^U)  >  are  zero.  We  also  show  that  the  only  non-zero  elements  of 
<Gmf(.t)G7iU)>  are  those  having  indices  which  come  in  pairs,  i.e.,  if  any  index  is  unique  and 
not  matched  by  some  other  index,  then  <Gmfi(t)Gyi(t)>  -  0.  Another  result  of  these  sym¬ 
metry  considerations  is  that  one  need  only  do  the  integration  over  the  sphere  for  two  adjacent 
octants,  saving  a  factor  of  four  in  computation  time. 

We  carry  out  our  computation  as  follows:  first  we  perform  the  integration  over  fl.  Using 
a  product  Gaussian  quadrature  method  for  numerical  integration  on  the  sphere,60  we  choose  a 
set  of  initial  directions  fl  for  the  scaled  initial  angular  velocity,  restricting  the  integration  to  two 
adjacent  octants.  For  each  direction  fl  we  compute  the  time  histories  of  the  elements 
Ga<,(fr,«q,fl)  using  Eqs.  (5.6)  and  (5.7)  and  a  variable  order,  variable  time  step  numerical 
integrator61  for  the  single  selected  uq,  adding  together  the  different  time  histories  corresponding 
to  different  fls  with  their  proper  weights,60  averaging  them  together  as  we  proceed  to  form  a 
single  set  of  time  histories  <?M(fr, «o)  for  the  infrared  case  and  Ga0 (f  r , «o)  G>g  (f  r , uq)  for  the 
Raman  case,  i.e.,  we  evaluate 

<?«,(*',« o)  -  / <?..(* t,uo,a)d(i  (5.20) 

and 

Ga(l(fr,uo)<?ra(fr.«o)  -  fGmfi((t,u0,{l)Gy>((t,u0.n)dn  .  (5.21) 

To  compute  the  final  <Gaa(t)>  and  < Gae(t)Grt(t)>  we  now  integrate  the  stored 
<G„a(ft,uo)>  and  <Ga0({f,uo)Gr4(f  r,uo)>  over  (  as  shown  in  Eq.  (5.19). 

These  free  rotor  results  may  also  be  extended49*56  to  treat  various  rotational  diffusion 
models. 

In  Fig.  3  we  present  <Gon(/)>  for  a  -  x,y,z,  along  with  their  Fourier  transforms,  for 
HjO.  All  relevant  parameters  for  the  molecule  are  calculated  ab  initio,  the  moment  of  inertia 
tensor  being  derived  from  the  ab  initio  equilibrium  geometry.  The  fact  that  G„,  does  not  tend 
to  zero  for  large  t  indicates  that  a  water  molecule  can  rotate  for  long  periods  about  its  smallest 
moment  of  inertia,  which  is  then  reflected  in  the  delta  function  Q  branch  seen  in  its  Fourier 
transform,  also  shown  in  Fig.  3. 

VI.  INFRARED  SPECTRA 

We  now  derive  the  appropriate  equations  for  ab  initio  computation  of  infrared  spectral 
band  contours.  While  nothing  in  this  section  is  basically  new  (for  example  see  the  work  of 
Gordon,16*11  of  Bratos,  Guissani,  and  Leicknam,49*56  and  our  earlier  papers,1*2  ),  we  derive 
here  in  one  place  the  various  equations  needed  for  the  different  approaches  discussed  in  this 
paper.  We  first  obtain  equations  appropriate  for  a  simultaneous  full  quantum  force  classical  tra¬ 
jectory  (QFCT)  spectral  computation,  and  then  derive  equations  for  the  separation  of  vibration 
and  rotation  and  the  treatment  of  rotation  in  a  fully  classical  rigid  rotor  approximation,  for  both 
simultaneous  and  sequential  approaches.  Finally  we  show  the  equations  for  the  sequential  rigid 
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rotor,  normal  mode  approximation. 

We  start  with  the  following  linear  response  theory  equations  from  Gordon11618  giving 
the  equilibrium  infrared  spectrum  as 

«<«)  -  4ir— ,(w)  •  C6.1> 

Wfcn 

m 

I(u)  —  (2ir)~xJ  dt  e~M <£(0)*/i(/)>  ,  (6.2) 

in  which  a(u)  is  the  absorption  cross  section  as  a  function  of  frequency  u  (not  to  be  confused 
with  the  angular  velocity  vector  «*),  P“(kg  T)~l  in  which  kg  is  Boltzmann’s  constant  and  T  the 
temperature,  1l—h/2ir  in  which  h  is  Planck's  constant,  c  is  the  speed  of  light,  n  is  the  index  of 
refraction  of  the  medium,  Hu)  is  defined  as  the  absorption  lineshape  and  is  evaluated  for  an 
isotropic  medium,  and  <£«))■£(/)>  is  the  ensemble  average  of  the  space  fixed  dipole 
moment  time  autocorrelation  function.  We  now  use  tildes  over  the  dipole  moment  vectors  fi, 
to  differentiate  these  space  fixed  dipole  moments  from  the  body  fixed  dipole  moments  to  fol¬ 
low,  which  we  will  write  without  the  tilde,  as  p.  Denoting  the  Fourier  transform  by  F,  we  can 
write  Hu)  as 

Hu)  -  (2w)-,/‘(<m(0)  mO)>}  ,  (6  3) 

where  we  define  the  Fourier  transform  as 

FfO)  -  J  e~,m,f(t)dt .  (6.4) 

We  can  also  express  Eq.  (6.2)  in  terms  of  the  power  spectrum1 


Hu)  —  (2w)~l<lim 


where  an  average  is  assumed  over  a  proper  ensemble  of  trajectories.  Given  an  ensemble  of 
time  histories  of  £(r)from  simultaneous  QFCT  calculations,  we  can  use  Eqs.  (6.1)  and  (6.2)  or 
(6.S)  to  compute  the  infrared  spectrum.1 

We  now  take  the  rigid  rotor,  vibrator  approximation  that  the  dipole  moment  can  be  writ¬ 
ten  as  £(/)  -  D(/)#t(r),  where  jkO)  is  the  space  fixed  dipole  moment,  DO)  is  the  rotation 
matrix  transforming  a  vector  relative  to  the  body  fixed  frame  to  a  vector  relative  to  space  fixed 
axes,  and  fi(/)  is  the  dipole  moment  in  the  body  frame.  Expressing  this  in  terms  of  com¬ 
ponents,  we  have  (again  using  the  summation  convention  that  repeated  indices  on  different 
symbols  are  to  be  summed  over) 

iktO)  —  DiaO)fiaU) ,  (6.6) 

Hu)  -  (2ir)-'F(<Dla(0)Ma(0)D,gO)tigO)>}  .  (6.7) 

We  note  that  because  of  the  linearity  of  both  the  Fourier  transform  and  the  ensemble  average, 
we  are  free  to  take  the  implied  summations  in  Eq.  (6.7)  wherever  we  wish. 

We  take  the  rigid  rotor,  vibrator  approximation  that  vibration  and  rotation  are  uncoupled 
and  uncorrelated,  and  thus  we  may  take  their  ensemble  averages  separately 

Hu)  -  (2ir)~,Fl<Dla(0)DigO)><fi„(0)ixgO)>)  .  (6.8) 

In  order  to  take  advantage  of  the  rotational  isotropy  of  space  as  discussed  in  Section  V 

above,  we  write  D(r)  as 

D(r)  -  D(0)G(r)  .  (6.9) 

We  now  have 


Hu)  -  (2w)_l  F[<Dia(0)Dly(0) GygO)>  <na(0)fig0) > I 


(6.10) 
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Because  D  is  orthogonal,  DlaDly  -  5ay,  and  we  have 

/(«)  -  (2»)-t  F { <  GtfU)  >  </*« (0)nf(f)  > )  .  (6.11) 

As  we  show  in  Appendix  B,  the  off-diagonal  components  of  <  Gag(.t)>  are  zero.  Accordingly, 
we  have  (making  the  summation  explicit  now) 

/(«)-(2w)-'  I  /,{<G„(f)><M.(0)M.(/)>}.  (6.12) 

This  can  also  be  written  in  the  following  two  ways:62 


/(«)-( 2tr)-'  £  /,{<Goa(r)>}*/’(<Mo(0)Ma(f)>) 

m—Xjra 


and 


/(«)  —  (2sr)->  £  F{<CM(r)>)*< 

•-xja 


V.0> 


>  , 


(6.13) 

(6.14) 


where  *  denotes  convolution  defined  as 


m 

f*g  («)  -  (2ir)-‘ J7(w  -  c-’)r(«')rf«' 


(6.15) 


Eqs.  (6.1)  and  (6.13)  or  (6.14)  am  appropriate  for  a  simultaneous  QFCT  rigid  rotor,  vibra¬ 
tor  approach  in  which  we  compute  the  classical  vibrational  trajectories  n„(t)  in  body  fixed  coor¬ 
dinates  and  then  convolve  with  the  rigid  rotor  rotational  band  shapes  f'{<(7„a(r)>).  In  order 
to  carry  out  a  sequential  rigid  rotor,  normal  mode  analysis,  in  Appendix  C  we  evaluate  the 
vibrational  correlation  function  F{<M<»(0)Ma(f)>}  in  terms  of  normal  coordinate  derivatives. 
Using  those  results,  we  obtain 


/(«) 


-  I  F{<Gtta(/)>}* 


(tiysu)  +  it^)2 

j~\ 


X  Huj  —  a») 

2m j  1  -  txp(-pXm)  ’ 


(6.16) 


where  is  the  equilibrium  value  of  fia,  8(m)  is  the  Dirac  delta  function,  qt  is  a  normal  coordi¬ 
nate,  and  wj  is  the  frequency  of  the  jth  normal  mode.  The  factor  to  the  left  of  the  *  in  the 
convolution  is  responsible  for  the  bandshapes,  while  the  factor  in  brackets  to  the  right  specifies 
the  frequencies  and  intensities  of  the  pure  rotational  and  vibrational  bands. 

It  is  necessary  to  apply  a  detailed  balance  quantum  correction1  to  the  classical  rotational 
results  for  all  the  approaches  given  above.  The  infrared  rotational  and  vibrational-rotational 
bands  for  equilibrium  systems  thus  are  corrected  by  multiplying  the  intensities  by 


exp(07r(Aw)/2] , 


(6.17) 


in  which  ft  —  ( k§T)~\  kg  is  Boltzmann’s  constant  and  Am  is  the  displacement  from  the  rota¬ 
tionless  biutd  center  (Aw  —  u  for  pure  rotational).  No  quantum  corrections  are  necessary  for 
the  vibrational  intensities  when  they  are  calculated  quantum  mechanically  as  in  Eq.  (6.16).  If 
the  vibrational  motion  is  treated  as  a  classical  trajectory,  as  in  Eqs.  (6.3),  (6.5)  or  (6.14),  then 
the  vibrational  rotational  band  intensities  can  be  multiplied  by  the  harmonic  quantum  correction 
factor1 


B*m 

1  -  exp(-£7T«)  ’ 

part  of  which  cancels  out  a  similar  factor  in  Eq.  (6.1). 


(6.18) 
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VII.  RAMAN  SPECTRA 

We  proceed  in  the  Raman  case  as  we  did  for  the  infrared.  First  we  derive  equations 
appropriate  for  a  simultaneous  QFCT  calculation,  then  from  this  we  obtain  equations  appropriate 
for  a  separation  of  rotational  from  vibrational  motion  and  finally  we  derive  the  equations  for  the 
more  usual  rigid  rotor,  normal  mode  sequential  calculation  following  the  lines  described  by  Bra* 
tos,  Guissani,  and  Leicknam.49**6  We  begin  with  the  linear  response  equations  for  Raman  spec- 
tn*  14. 16. 18 

-^L_  |  _  (2sr)-*jr  dt  exp(-/»r)|<Tr[?i“(0)Pto(r)j>  (7.1) 

(*'**  )  "  {2v)~XS  *  exp(-/W) <Tr[f- (0)f  (/)]  >  ,  (7.2) 

in  which  Eqs.  (7.1)  and  (7.2)  are  the  differential  cross  sections  for  scattering  into  angular  fre¬ 
quency  range  dot  and  solid  angle  range  d(l  for  the  isotropic  and  anisotropic  Raman  spectra, 
respectively,  weighted  by  in  which  2irit/  is  the  wavelength  of  the  incident  radiation,  and  by 
It,,  in  which  2ifXt  is  the  wavelength  of  the  scattered  radiation.63  P  is  the  space  fixed  Cartesian 
polarizability  tensor,  P1"  -  P**  1  is  its  isotropic  part,  where  F”  -  }(TrP)  and  1  is  the  identity 
tensor,  P"  -P  - P*  is  its  anisotropic  part,  Tr  indicates  trace,  and  it  is  understood  that  the 
correlation  functions  in  Eqs.  (7.1)  and  (7.2)  are  averaged  over  the  desired  ensemble.  We  now 
use  tildes  over  the  polarizability  tensors  P  to  differentiate  these  space  fixed  polarizabilities  from 
the  body  fixed  polarizabilities  to  follow,  which  we  will  write  without  the  tilde,  as  P.  We  can 
also  express  Eqs.  (7.1)  and  (7.2)  in  terms  of  the  power  spectrum2 

(*'*■’  -£fc  )„  (7J> 

H  L  ■  >  •  (7  4) 

where  averaging  over  the  proper  ensemble  of  trajectories  is  assumed.  Given  an  ensemble  of 
P(r)  time  histories  from  simultaneous  QFCT  calculations,  we  can  use  Eqs.  (7.1)  and  (7.2)  or 
(7.3)  and  (7.4)  to  compute  the  Raman  spectrum. 

For  both  the  isotropic  and  anisotropic  cases  we  wish  to  evaluate  time  correlation  functions 
of  the  form  <Tr(P(0)P(f)]>.  Employing  the  summation  convention  once  again,  we  write 


<Trt* (0)f(r)]>  -  <fl*(0 )**,(/)>  .  (7.5) 

As  in  the  infrared  case,  we  take  the  rigid  rotor,  vibrator  approximation  and  write 

P*0)  -  Dlait)Pa$(i)Drf(t)  ,  (7.6) 

where  P  is  the  space  fixed  and  P  is  the  body  fixed  polarizability  tensor.  Then 

<Tr(P(0)?(f)]>  -  <Dla(O)Pa0(O)Df}(O)Dkr(t)PYi(t)Diri'(t)>.  (7.7) 

Separating  vibration  and  rotation  as  before,  we  have 

<Trtf(0)*(/)]>  -  <Dla(0)Dfi(0)Dky(t)D;lHt)><Pae(0)Prt(t)>  .  (7.8) 


As  in  the  infrared  case,  we  may  write  D(r)  -  D(O)GO),  where  G(/)  is  the  rotational  Green's 
Function.  Making  use  of  the  fact  that,  because  rotation  matrices  are  orthogonal, 
D$jHt)  -  DJt(t),  we  write 

<Trl*(0)*(/)]>  -  < Z)/, (0) Dji {0)Dktl{0)GHV(t)D„X0)G,^t)>< P„g (0)PyS(t)> 

-  <  Dta  (0) Drf  (0)  (0) A-<0)  >  <  G„  (f )  G*  (f )  >  <  Pa„  (0)  Pyh(t)>  .  (7.9) 

Since  D~’(0)D(0)  -  1,  the  unit  matrix,  Drf(0)DiM(0)  -  6^,  and  because  D  is  orthogonal, 
Dim(0)Dj0)  —  5n,„  and  thus  we  have 
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R 


P 


P 


<TrtP(0)P(r)]>-  <  G$y (t) Ga6(t) >  < Pa0(O)PYt(t)  >  . 

We  now  make  use  of  the  relation  Pa,-P"°  *«  +  P3  to  write  Eqs.  (7.1)  and  (7.2)  as 


<P<r 

dud  Cl 

d2g- 
dud  Cl 


(6ir)-'Fl  < GfrU) Gat(f) > < F" (O)8a0F" (/)8y*>  ) 
(2ir)-'F{  <G0y(t)Ga,(t)><P%(O)/>fl(t)>  }  . 


(7.10) 

(7.11) 

(7.12) 


We  can  also  write  these  as 
<P<r 


dud  Cl 
d2<r 


(6 »)-'/•{  <  <?*(/)  Grt(/)>)  */'{</>i,B(0)8a(|/*<*#(f)8r»>  )  . 


(*!*/ 

jx,*/  J  -  (2 »)-'/•{  <G^(r)Ga4(r)>)  */•{</»& (0)P2(f)>  )  • 


and62 


(7.13) 

(7.14) 


(*'*•’  ~£dn  L  "  <o»U)G*M> I  * 

<  limyj-  .  (7.15) 

(*.*?  I  -  (2#)-'F[  <C#,(l)G.J(/)>l  • 


where  the  t  is  used  here  to  indicate  complex  conjugation. 

Eqs.  (7.13)  Mid  (7.14)  or  (7.1S)  and  (7.16)  are  appropriate  for  a  simultaneous  QFCT  tra¬ 
jectory  calculation  in  which  the  vibrational  part  is  computed  from  P(r)  along  trajectories  Tor 
rotationless  molecules,  and  the  band  contours  are  computed  from  F[<Gpy(t)GaiU)>\.  For 
our  simple  normal  mode  analysis,  in  Appendix  C  we  evaluate  the  vibrational  correlation  func¬ 
tion  /,{</>.,(0)Py,(/)>}  in  terms  of  normal  coordinate  polarizability  derivatives.  Using  those 
results,  we  obtain 


I* 


d2<r 
dud  Cl 


a.)-'  *  iUpf  JL  T-i 

,Ti|  »«/  2“/  1  - 


&(ujj-  (l>) 


exp(-pTlu) 


(7.17) 


I nrrtrtM  *  tp-  *p  -5-  8(“/-"> 

I  r>  H  bqt  2uj  1  - 


exp(-0ftw) 


(7.18) 


where  the  superscript  0  indicates  the  value  at  the  equilibrium  position,  qf  is  the  jth  normal 
coordinate,  and  uj  is  the  frequency  of  the  jth  normal  mode. 

Quantum  corrections  to  Raman  spectra2  can  be  made  exactly  as  for  IR  spectra  (see  the 
discussion  at  the  end  of  Section  V!).  Eqs.  (7.1)  to  (7.4)  and  Eqs.  (7.1S)  and  (7.16)  can  be 
quantum  corrected2  for  rotation  and  harmonic  vibration  by  Eqs.  (6.17)  and  (6.18),  while  Eqs. 
(7.17)  and  (7.18)  require  quantum  correction  only  of  the  rotational  motion,  using  Eq.  (6.17). 


VIII.  SEQUENTIAL  INFRARED  AND  RAMAN  APPROACH  FOR  HjO 

We  now  compute  the  infrared  band  contours  for  dilute  water  vapor  as  an  example  of  an 
ab  initio  sequential  approach  using  the  data  from  Tables  I  to  III  from  Sections  III  and  IV  and 
Eqs.  (6.1),  (6.16),  and  (6.17)  along  with  the  <Gaa(/)>  elements  evaluated  as  shown  in  Sec¬ 
tion  V.  Rotation  is  computed  classically  from  the  quantum  derived  equilibrium  moments  of 
inertia  and  vibration  is  treated  quantum  mechanically  in  the  normal  mode  (harmonic)  approxi¬ 
mation  to  produce  the  infrared  spectral  band  contours  for  the  water  molecule  in  the  lowest  rigid 
rotor,  normal  mode  approximation.  Figure  4  shows  the  infrared  spectral  band  contours.  The 
upper  trace  is  generated  by  broadening  the  water  lines  from  the  Air  Force  Geophysics  Labora¬ 
tory  (AFGL)  tape64’65  until  they  merge  into  band  contours.  The  AFGL  data  is  derived  from 
experimental  measurements  with  considerable  help  from  theoretical  calculations.  The  lower 
trace  is  our  simple  ab  initio  rigid  rotor,  normal  mode  computation  of  the  infrared  band  contours. 

The  presence  of  a  Q  branch  in  the  infrared  spectrum  depends  upon  the  value  of 
f(<C„>)  at  »  — 0.  The  value  of  a  Fourier  transform  at  w  -  0  depends  simply  on  the 
integral  over  all  time  of  the  untransformed  function.  If  this  area  is  infinite,  then  the  Fourier 
transform  has  a  delta  function  at  at  -  0.  A  delta  function  in  F  <G„>  becomes  a  Q  branch 
when  put  in  Eq.  (6.16)  if  there  is  a  nonzero  dipole  derivative  in  the  a  direction.  The  area 
under  <Gaa(t)>  is  infinite  if 

lim<G««(/)>  a*  0  .  (8.1) 

f— •• 

<Goa(f)>  can  be  viewed  as  the  ensemble  averaged  autocorrelation  function  of  a  unit  vector 
fixed  in  the  molecule  in  the  direction  of  the  a  principle  axis.  If  rotation  about  one  of  the  prin¬ 
ciple  axes  is  very  stable,  then  condition  (8.1)  can  hold,  and  there  will  be  an  infrared  Q  branch 
associated  with  each  normal  mode  that  has  a  dipole  derivative  component  in  the  direction  of 
that  principle  axis.  In  general,  rotation  about  intermediate  moments  of  inertia  is  very 
unstable,39  so  that  one  would  not  expect49  any  Q  branches  associated  with  them.  Rotation 
about  the  large  and  the  small  moments  of  inertia  is  generally  more  stable.  Asymptotic  formulae 
for  <(rM(/)>  are  discussed  by  Guissani,  Leicknam,  and  Bratos.49  For  the  case  of  H2O,  rota¬ 
tion  about  the  smallest  moment  of  inertia,  the  y  axis  as  shown  in  Fig.  1,  satisfies  condition 
(8.1),  as  shown  in  Fig.  3.  We  obtain  a  limiting  value  for  that  axis  of  0.344  and  thus  the 
Fourier  transform  has  a  delta  function  of  area  2vx  0.344.  Since,  as  shown  in  Table  II|,  the  v} 
vibration  has  a  dipole  derivative  component  in  the  y  direction,  we  obtain  a  Q  branch  centered 
at  vy  As  can  be  seen  from  comparing  with  the  experimental  spectra  in  Fig.  4,  this  rigid  rotor 
central  Q  branch  has  an  intensity  that  is  too  large.  We  surmise  that  various  effects  present  in 
real  systems,  but  not  in  our  ideal  rigid  rotor,  serve  to  alter  the  long  term  behavior  of 
<  G„(t)>  and  thus  the  shape  of  the  Q  branch  in  the  real  spectrum.  For  instance,  the  coupling 
between  vibration  and  rotation  or  collisions  could  serve  to  disturb  the  ordered  rotational  effects 
which  lead  to  condition  (8.1).  This  would  allow  <  G„(t)>  eventually  to  die  off,  resulting  in  a 
less  intense  Q  branch  more  in  accord  with  the  experimental  spectrum.  The  v\  symmetric 
stretch  vibrational  rotational  band  is  very  small  in  comparison  with  the  vy  band,  and  is  hardly 
seen.  The  v2  symmetric  bend  band  has  no  Q  branch  because  its  non-zero  dipole  moment 
derivative  is  along  the  intermediate  moment  of  inertia  axis  and  rotation  about  this  axis  is 
unstable.  Hie  difference  between  the  theoretical  and  experimental  ratios  of  the  P  and  R  type 
branches  in  the  vj  band  is  at  the  moment  unexplained.  For  a  different  approach,  in  which  an  ab 
initio  dipole  moment  function  is  used  to  compute  the  strengths  of  vibration  •  rotation  lines  see 
the  work  of  Carney.66 

Figure  5  shows  the  isotropic  and  anisotropic  rigid  rotor,  normal  mode  Raman  spectral 
band  contours  of  dilute  water  vapor  as  given  by  our  ab  initio  calculations,  using  the  data  from 
Tables  Mil  and  Eqs.  (7.17),  (7.18),  and  (6.17)  with  the  elements  of  <Gpy(t)Gai(t)>  com¬ 
puted  as  shown  in  Section  V.  Note  that  in  our  computed  Raman  spectra,  the  OH  stretching 
region  has  strong  contributions  from  both  the  totally  symmetric  stretch  v\  and  the  asymmetric, 
stretch  vy ,  as  is  seen  experimentally,46  while  the  infrared  stretching  region  is  dominated  by  the 
asymmetric  stretch  vy.  The  results  for  the  isotropic  case  from  our  calculations  are  actually  delta 
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Finn*  4.  Comparison  of  measured  and  computed  infrared  spectral  band  contours  for  H20  va¬ 
por  at  296  K.  The  ’experimental"  spectrum  is  derived  by  broadening  the  spectral  lines  from  the 
Air  Force  Geophysics  Laboratory  tape,  and  is  plotted  with  its  baseline  raised  by 
15x10-*  m2  molecule''  for  clarity.  The  ab  initio  sequential  theoretical  spectrum  is  computed  in 
the  rigid  rotor,  normal  mode  approximation  using  quantum  mechanical  normal  mode  vibration 
with  a  calculation  (from  the  quantum  derived  equilibrium  moments  of  inertia)  of  the 

band  shapes  arising  from  rotational  motion.  The  left  and  tallest  peak  is  pure  rotation  the  mid- 
dtepeaksat  approximately  1500  to  2000  cm"1  are  the  v2  bending  mode  and  the  peaks  on  the 
right  ^approximately  3500  to  4500  cm~'  are  the  OH  stretching,  the  «m  symmetric  stretch  being 
so  weak  as  to  be  hardly  seen  and  the  vj  asymmetric  stretch  *hus  dominating  the  stretching  re¬ 
gion.  The  sharp  peak  in  the  stretching  region  of  the  theoretical  spectrum  is  the  i>i  Q  branch, 
which  can  be  seen  to  be  much  broader  and  weaker  in  the  experimental  spectrum,  indicating 
that  the  real  long  term  vectorial  correlation  dies  oflT  more  quickly  than  the  rigid  rotor  one.  Note 
that  the  theoretical  vibrational  rotational  peaks  occur  at  higher  frequencies  than  the  experimen¬ 
tal  ones  because  anharmonicity  was  not  included  in  this  normal  mode  calculation. 
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Figure  S.  Ab  initio  sequential  isotropic  and  anisotropic  Raman  spectral  band  contours  Tor  HP 
vapor.  The  spectra  are  computed  in  the  rigid  rotor,  normal  mode  approximation  using  quantum 
mechanical  normal  mode  vibration  with  a  classical  rotational  calculation  of  the  band  contours 
from  the  quantum  equilibrium  moments  of  inertia.  The  temperature  is  2%  K.  The  overall  band 
contour  in  the  Raman  OH  stretch  region,  3500  to  4500  cm~\  has  significant  contributions  from 
both  the  symmetric  stretch  v\  and  the  asymmetric  stretch  vj,  in  contrast  to  the  infrared  OH 
stretch  region,  which  is  dominated  by  v3.  The  sharp  Q  branch  peak  in  the  Raman  stretching  re¬ 
gion  is  due  to  vi,  while  that  for  the  infrared  region  is  due  to  i/3.  Due  to  space  limitations,  the 
heights  of  the  anisotropic  Q  branches  for  pure  rotation  and  for  v\  are  cut  ofT  at  the  upper  mar¬ 
gin  of  the  figure.  The  anisotropic  spectrum  is  plotted  with  a  vertical  offset  of 
1.5xl0-76  m2ster~\  in  which  ster  means  steradians.  Note  that  the  experimental  vibrational  rota¬ 
tional  tends  will  be  somewhat  lower  in  frequency  than  the  normal  mode  frequencies  shown 
here  due  to  vibrational  anharmonicity. 
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functions  at  0,  kj,  and  v2.  There  is  nothing  for  in  the  isotropic  spectrum  because  the 
corresponding  isotropic  polarizability  derivatives  along  the  diagonal  in  Table  III  are  zero.  We 
have  broadened  the  delta  functions  so  as  to  be  able  to  present  them  graphically.  Although 
experimental  Raman  line  spectra  of  water  vapor  have  been  measured,46-48*67’68  the  available 
experimental  band  contours  67  are  of  low  signal  to  noise  ratio  and  make  accurate  comparison 
with  our  calculations  difficult.  Water  vapor  pressurized  by  an  inert  gas  could  provide  lines 
appropriately  broadened  and  merged  into  a  band  contour  for  an  experimental  comparison  with 
this  calculation.  A  comparison  with  quantum  theory  could  be  made  by  broadening  the  line 
spectrum  from  an  accurate  quantum  calculation.47’48 

Several  points  can,  however,  be  deduced  from  the  available  Raman  band  contour  and  line 
spectral  data. 46-48*67’ 68  As  is  also  seen  experimentally,48'67  our  pure  rotational  anisotropic  band 
at  approximately  0  to  500  cm~x  is  indeed  much  weaker  than  the  corresponding  pure  rotational 
band2  for  N2.  Our  computed  equilibrium  polarizability  matrix  elements  given  in  Table  II,  which 
produce  the  rotational  spectrum,  compare  quite  well  with  those  deduced  experimentally  by 
Murphy,48  when  his  axis  convention  is  converted  to  ours42  by  x— >z,  y  —  x,  z—  y.  The  com¬ 
puted  ?2  bending  band  at  approximately  1500  to  2000  cm~x  is  indeed  much  weaker  than  the  OH 
stretch  region,  as  is  found  experimentally,47-67  and,  as  well  as  we  can  judge  from  the  published 
figures,  the  computed  *2  band  shape  at  least  qualitatively  resembles  that  of  the  observed  line 
spectrum,47  except  that  our  i>2  Q  branch  appears  somewhat  sharper  than  the  experimental  one. 
The  computed  OH  stretch  region  (the  3500  to  4500  cm~x  band)  also  resembles  the  available 
survey  band  spectrum  of  Penney  and  Lapp67  and,  as  well  as  we  can  judge  from  the  published 
figures,  also  resembles  the  line  spectrum  measured  by  Murphy.46  In  Murphy's  line  spectrum46 
we  see  the  strong  v\  Q  branch  as  in  our  spectrum  (although  the  experimental  Q  branch  is 
broader),  as  well  as  the  major  contribution  from  v3  with  no  Q  branch  as  we  also  observe  in  our 
calculation. 

Note  that  the  computed  infrared  and  Raman  vibrational  rotational  bands  are  at  higher  fre¬ 
quencies  than  the  experimental  ones  because  only  the  harmonic  forces  were  used  in  this  simple 
normal  mode  calculation,  and  the  anharmonicity  which  we  did  not  compute  lowers  the  real  fre¬ 
quencies. 

IX.  SIMULTANEOUS  QUANTUM  FORCE  CLASSICAL  TRAJECTORY  APPROACH 

In  addition  to  the  sequential  approach  illustrated  in  the  H20  example  presented  in  Section 
VIII  above,  a  simultaneous  quantum  force  classical  trajectory  (QFCT)  approach  can  also  be 
taken  which  is  likely  t'.  be  more  practical  for  many  atom  systems.  Given  the  force  F,  on  each 
atom  i  at  position  r,  from  the  quantum  gradient  technique  described  in  Section  III,  we  can  in 
principle  integrate  forward  one  step  in  time  the  classical  equations  of  motion 

d\ 

F,  -  "  1 . N  (9.1) 

for  the  N  atoms.  (For  nearly  harmonic  systems  such  as  isolated  molecules,  the  trigonometric 
polynomial  fitting  routine  of  Gautschi69  may  be  of  interest.)  Then,  given  the  new  set  of  posi¬ 
tions,  we  can  derive  from  the  quantum  part  of  our  calculation  new  values  of  the  F,,  and  of 
whatever  quantities  such  as  p  (for  the  infrared  spectrum),  and  P  (for  the  Raman  spectrum)  we 
wish  in  addition. 

From  an  ensemble  of  the  time  histories  #*(/)  and  P(r),  we  can  derive  the  infrared1  and 
Raman2  spectra,  using  Eqs.  (6.1)  to  (6.5)  and  (7.1)  to  (7.4),  in  either  the  correlation  function 
or  power  spectrum  forms.  The  appropriate  simple  quantum  corrections  shown  in  Eqs.  (6.17) 
and  (6.18)  are  discussed  in  other  papers.1 -2’ 14  As  illustrated  elsewhere,  these  techniques  can  be 
extended  to  calculate  the  infrared  and  Raman  spectra  of  non-equilibrium  and  time  dependent 
systems,14’70  as  well  as  the  electronic  absorption  spectra  of  equilibrium  and  non-equilibrium 
systems.71 

Because  the  computational  demands  of  a  many  atom  simultaneous  approach  are  high,  care¬ 
ful  thought  is  warranted  as  to  how  to  efficiently  average  over  the  ensemble  of  interest.  For  an 
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equilibrium  ensemble,  we  can  start  a  series  of  runs  by  approximating  the  initial  phase  space  dis¬ 
tribution  in  the  following  manner.  If  the  system  in  question  were  completely  harmonic,  on  the 
average  its  total  internal  energy  would  be  equally  divided  into  kinetic  and  potential  energy.  If 
we  start  our  atoms  at  the  equilibrium  positions,  but  give  them  velocities  chosen  randomly  from 
a  Boltzmann  distribution  corresponding  to  twice  the  temperature,  we  would  give  the  molecule 
the  correct  distribution  of  total  energy  for  the  true  temperature  if  it  were  truly  harmonic.  One 
could  remove  angular  momentum  by  adding  on  a  rotational  term  to  the  velocities  correspond¬ 
ing  to  a  coordinate  system  rotating  in  the  opposite  direction  to  the  initial  angular  velocity,  and 
then  rechoose  a  rotational  angular  momentum  for  the  system  corresponding  to  the  real  tem¬ 
perature  of  interest.  For  an  anharmonic  system  this  is  only  an  approximation;  in  addition, 
starting  each  run  at  the  equilibrium  position  for  even  a  multidimensional  harmonic  oscillator 
may  not  sample  all  of  phase  space.  One  therefore  would  use  this  method  only  for  the  first 
equilibration  run.  For  each  run  thereafter,  one  would  use  the  final  positions  in  the  previous 
run  as  initial  positions,  and  give  a  new  set  of  random  initial  velocities  chosen  from  the 
Boltzmann  distribution  corresponding  to  the  true  temperature.  After  a  number  of  such  runs, 
the  system  would  be  equilibrated  at  the  desired  temperature,  and  one  could  begin  to  take  data. 

An  advantage  of  the  simultaneous  QFCT  ab  initio  molecular  dynamic  approach  described 
above  is  its  simplicity.  No  separation  of  translational,  rotational,  and  vibrational  motion  and  no 
normal  mode  approximations  are  needed.  The  computation  may  equally  well  be  carried  out  for 
gases  of  individual  molecules,  clusters,  liquids,  solutions,  or  solids.  Many  other  measurables 
besides  infrared,  Raman,  and  electronic  spectra  may  also  be  computed  in  a  parallel  QFCT 
manner,  by  quantum  mechanically  calculating  parameters  along  the  trajectory  and  then  per¬ 
forming  the  appropriate  averages  over  an  ensemble  of  trajectories.  For  example,  thermo¬ 
dynamic  variables  such  as  energy,  heat  capacity,  free  energy,  and  entropy  can  be  computed  in  a 
QFCT  manner  from  quantum  evaluation  of  the  potential  energy  seen  by  the  nuclei  and  classical 
evaluation  of  the  nuclear  kinetic  energy  with  averaging  over  proper  ensembles  of  trajectories.19 
In  addition,  by  keeping  track  of  velocities  along  the  classical  paths  and  power  spectral  transfor¬ 
mation  to  velocity  spectra,  the  quantum  corrections  to  these  thermodynamic  variables  may  also 
be  computed.19 

In  cases  in  which  the  above  approach  is  too  demanding  in  terms  of  available  computer 
power,  again,  as  in  the  sequential  cast,  a  rigid  rotor  vibrator  approximation  can  be  made  and  the 
effect  of  the  rotational  motions  handled  entirely  classically  once  the  quantum  equilibrium 
geometry  is  known.  For  example,  to  compute  infrared  spectra  in  this  way  one  would  use  Eqs. 
(6.1),  (6.12)  to  (6.15),  (6.17)  and  (6.18)  and  to  compute  Raman  spectra  Eqs.  (7.11)  to  (7.16), 
(6.17)  and  (6.18),  choosing  either  the  convolution  or  power  spectral  methods. 

X.  DISCUSSION  AND  CONCLUSION 

We  have  shown  that  a  judicious  combination  of  quantum  mechanics  and  classical  mechan¬ 
ics  can  be  used  to  compute  ab  initio  infrared  and  Raman  spectral  band  contours.  Two  different 
approaches  are  discussed.  In  the  first,  one  sequentially  solves  the  electronic  and  then  the 
nuclear  parts  of  the  Bom-Oppenheimer  approximation,  and  in  the  second  they  are  solved  simul¬ 
taneously.  We  have  presented  an  example  at  the  simplest  rigid  rotor,  normal  mode  level  for  the 
sequential  approach.  The  accuracy  of  this  simplest  approach  is  not  exceptional,  but  it  could  cer¬ 
tainly  be  considerably  improved  by  bringing  in  effects  beyond  the  rigid  rotor,  normal  mode 
level  such  as  centrifugal  distortion,  vibration  rotation  interaction  (for  example  Coriolis  effects) 
and  anharmonic  effects.  In  addition,  the  effects  of  terms  beyond  the  linear  expansion  for  the 
dipole  moment  and  polarizability  tensor  could  be  considered. 

In  developing  the  methodology  for  the  sequential  rigid  rotor,  normal  mode  approach,  we 
have  worked  out  an  improved,  simple,  accurate,  and  fast  technique  to  compute  pure  rotational 
and  vibrational  •  rotational  classical  rigid  rotor  infrared  and  Raman  band  contours  for  any 
molecule,  including  asymmetric  tops.  Comparison  with  the  pioneering  work  on  asymmetric 
rotors  by  Leicknam,  Guissani,  and  Bratos49*56  shows  agreement  at  intermediate  and  higher  dis¬ 
placements  from  the  band  center,  but  some  quantitative  disagreement  nearer  the  band  centers 


where  we  believe  we  were  able  to  achieve  higher  accuracy.  The  overall  qualitative  picture  is  the 
same,  but  we  achieve  it  with  a  simple  and  rapid  calculation  which  avoids  complexities  such  as 
reference  to  Euler  angles  and  special  functions. 

We  discuss  a  second  simultaneous  approach,  which  we  call  Quantum  Force  Classical  Tra¬ 
jectory  (QFCT),  which  is  more  appropriate  for  many  atom  systems  such  as  large  molecules, 
clusters,  and  liquids.  In  this  latter  approach,  which  closely  follows  the  spirit  of  our  earlier 
molecular  dynamics  approaches  to  infrared1  and  Raman2  spectra,  there  is  no  need  to  separate 
vibrational,  rotational,  and  translational  motion  or  to  make  normal  mode  approximations.  The 
generality  of  the  approach  is  such  that  it  can  in  principle  be  applied  to  gas,  liquid  or  solid  sys¬ 
tems,  composed  of  stable  or  unstable  species,  in  equilibrium  or  not.14  The  QFCT  method  could 
also  be  employed  to  compute  other  spectral  (for  example  electronic  absorption71 )  and  transport 
properties  derivable  from  correlation  functions,  to  deduce  thermodynamic  quantities,19  and  to 
study  chemical  reactions,  as  has  already  been  demonstrated  by  Warshel  and  Karplus12  and  by 
Leforestier.13  A  parallel  Monte  Carlo  approach  in  which  energies  are  computed  quantum 
mechanically  at  each  trial  step  and  other  measurables  computed  at  the  accepted  steps  might  be 
used  in  the  computation  of  parameters  which  do  not  depend  on  momenta. 
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APPENDIX  A 

In  the  classical  statistical  mechanics  of  a  rigid  body,  we  are  required  to  calculate  equili¬ 
brium  averages  by  integration  over  phase  space,  using  the  canonical  volume  element 
dT  —  dqidq1dq3dp]dp2dpi,  where  qa(<r  —  1,2,3)  are  coordinates  in  the  rotation  group  (eg., 
Euler  angles)  and  p„  are  the  canonically  conjugate  momenta.  In  practice,  it  is  often  easier  to 
use  the  volume  element  dV  dl\dl2dl2,  where  /„  (m  —  1,2,3)  are  the  components  of  the  angular 
momentum  in  the  body  frame,  related  to  the  angular  velocity  by  /M  -  (/M„  is  the  moment 

of  inertia  tensor  in  the  body  frame  and  we  again  in  this  appendix  use  the  summation  notation), 
and  dV  is  the  invariant  volume  element  in  the  rotation  group.  We  have  not  found  an  adequate 
reference  for  this  procedure,  so  we  justify  it  here. 

The  connection  between  q  and  w  is  given  by  the  rotation  matrix  Dia(q)  [see  Eq.  (5.7)] 

A.  -  I*  -  "Sto.D^q)  ,  (Al) 

so  that  the  angular  velocity  u  in  the  body  frame,  is  related  to  the  configurational  velocity  q  by 
”  i*gM«A g(?>  ^  W  •  (A2) 

(The  second  equality  defines  the  matrix  A.)  The  Lagrangian  is 

I*  ™  \ '  (A3) 


so  the  canonical  momenta  are 


/  ,  j 


t,.A 


t'tT  * 


The  canonical  phase  space  volume  element  is 


dT  -  dqldq2dq3(det  A)<Hldl2dlJ 


(A4) 
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-  dV  dl\dl2dl}  (AS) 

with  dV  —  (del  A)dqxdq1dqi,  which,  we  claim,  is  the  invariant  volume  element  in  the  rotation 
group. 

Define  the  differential  forms  X„  -  Avadq,T,  so  that  dV  -  x>X2Xj  (exterior  product).  The 
forms  X„  are  left  invariant: 

X,  -  A^dq-  -  dd°  -  2e0*aDifi(q)dD,a(q).  (A6) 

For  fixed  r, 

\,(rq)  -  \^Dlfi(rq)dDim(rq) 

-  i«^y<r)D^)D*(r)dPka(9)  (A7) 

■  ■  X,(q)  ,  (A8) 

where  we  use  the  orthogonality  of  D*ir).  Since  the  X,  are  left  invariant,  so  is  their  exterior 
product  dV. 

The  probability  distribution  in  phase  space  is  dTtxp(—pH),  where  H  is  the  Hamiltonian. 
We  have 

H  -  .  (A9) 

and 

dT  -  dVdlidljdij  -  det(/)</K</u,</w2</«j ,  (AlO) 

so  the  averaging  procedure  of  Eq.  (S.10)  has  been  justified. 


APPENDIX  B:  SYMMETRY  IN  THE  INERTIA  TENSOR 

We  can  take  advantage  of  symmetry  in  the  molecule’s  inertia  tensor  to  show  that  most  of 
the  elements  and  products  of  <G>  are  zero.  Consider  a  rotation  R  which  commutes  with  the 
inertia  tensor  I,  so  that  IR  —  RI  .  We  will  now  show  that  if  w(r,«°)  satisfies  the  equations  of 
motion,  then  so  does  R«(/,a»°).  We  have,  using  the  summation  notation, 

Arf-^(R«*)a  ”  LpRfa—jf-  (BD 

-  <»» 

■  —Rafily *  (B3) 


in  which  is  again  the  Levi-Civita  density.  Since  R  is  orthogonal,  R„,.Rav-  -  6,v  and 


RryRjy  -  &yy,  SO  we  C8H  write 

fa 0  jWa)|J  "  Raffi**'  p‘y'0  (B4) 

“  R a&R <tvR<tv'R ryR ry'  ^ v'y'pl p*** v  •  (B5) 

Using  the  definition  of  a  determinant  and  the  fact  that  R  is  orthogonal  one  can  show  that 
Ra$RirrRry'^¥,y‘0  ™  *a<yr  det  R  ■  ta<rT  ,  (B6) 

where  det  is  the  determinant.  This  leads  to 

f«a-j(R*»>*  "  ~RapRry*aa rfyMW*.W*-  <®7) 

-  -*aaTlTy{ Ra)y(Ra),  (B8) 


I 


which  is  just  Eq.  (S.S)  with  R»  in  place  of  •».  Because  m(/,Rm°)  and  R satisfy  the 
same  equations  of  motion  and  have  the  same  initial  conditions,  they  must  be  equal,  so  that 
«t(f,R«*°)  -  R«0,«°). 

We  will  now  show  that  if  the  set  D (t),m(t)  satisfies  the  equations  of  motion,  then  so  does 
D(/)R~l,R«*(f).  Putting  D(/)R-1  into  Eq.  (S.7),  we  have 


Ra0UfiDjyRTy 


-  «MT  (Rw)»(DR-,)<r  .  (B12) 

For  any  given  solution  of  Eqs.  (S.6),  the  solution  of  Eq.  (S.7)  is  unique,  for  given  initial  condi¬ 
tions.  When  we  are  dealing  with  G,  however,  the  initial  conditions  are  always  the  same: 

G(0,»°)  —  1.  By  definition,  G(r,R«°)  is  a  solution  to  Eq.  (S.7)  given  that  the  solution  to  Eqs. 
(S.6)  is  Rat0.  Similarly,  we  have  just  obtained  that  D(/)R_1  is  a  solution  subject  to  the  same 
solution  for  at.  The  initial  conditions  do  not  match  however,  but  since  D(/)R_I  is  a  solution, 
by  the  isotropy  of  space,  so  is  RO(/)R't,  which  matches  initial  conditions.  This  implies  that 

GO,Rat°)  -  RG(/,at°)R-'  .  (B13) 

For  the  infrared  case  49'50  we  can  use  Eq.  (B13)  to  show  that  <G(r)>  is  diagonal.  I  in 
the  body  frame  commutes  with  each  of  the  following  180*  rotations: 

[  1  0  01  [-10  0]  [-1  0  01 

-  0  -1  0  ;  <r,  -  0  1  0  ,  «rt  -  0  -1  0  .  (B14) 

l  o  o  -ij  ^  loo  -11  l  0  0  l| 

For  any  off  diagonal  element  of  G,  G„g  with  a^/3,  (voG®’“,)(,0  -  -Ga0  (no  summation),  so 
putting  this  into  Eq.  (B13) 

Ga0(t,i ram°)  -  —Ga0(t,m°)  .  (B15) 

For  a  Boltzmann  distribution  of  angular  velocities  equal  weight  is  given  to  m  and  «•„•»,  so 
<Ga0(t)>  -  -<<?„*(*)>  -  0,  a*/3  ,  (B16) 

and  we  conclude  that  <GO)>  is  diagonal. 

In  order  to  calculate  Raman  spectra  we  need  to  calculate  <GagGyh>.  If  some  one  of 
(a,/9,-y,8)  is  unique,  say  /9  (i.e.,  p*y,  P*a,  and  0*8),  then  application  of  <r0  as  above  will 
change  the  sign  of  Ga0Gyi  and  we  will  get  <  G„gGyi>  -  0.  Therefore  the  indices  must  come 
in  pairs,  and  the  nonzero  elements  are  of  the  form  <  G*g  >,  <  Ga0G0a>,  or  <  GaaGgg>. 

Because  with  the  application  of  each  of  the  rotations  in  Eq.  (B14)  one  can  generate  the 
entire  sphere  from  any  two  adjacent  octants,  and  one  need  only  do  the  averaging  over  the 
sphere  for  two  adjacent  octants. 

APPENDIX  C:  DERIVATION  OF  THE  VIBRATIONAL  CORRELATION  FUNCTIONS 
Given  two  observables  A  and  B,  we  wish  to  evaluate 

F\<A  ((»*(/)>}  (Cl) 

in  terms  of  normal  coordinates.  In  this  appendix  we  will  not  use  the  summation  notation.  For 
the  infrared  case  we  take  A  -  B  -  /*«»  and  for  the  Raman  case  A  -  Pa0  and  B  -  Pyt.  We 
approximate  A  and  B  in  a  Taylor’s  series  in  terms  of  normal  coordinates  qJy  keeping  only  terms 
to  first  order,  obtaining 

A  (0)  -  A0  +  A4  -  A0  +  Zp-Qj  (0)  (C2) 
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BU)  -  fl°  +  lfl(r)  -  fl°  +  (C3) 

j  9*j 

where  A9  and  denote  the  values  of  A  and  B  at  equilibrium,  and  the  derivatives  are 
evaluated  at  equilibrium. 

From  the  fluctuation-dissipation  theorem11  and  the  Kubo  formula,1*  we  have 

/*{<i4(0)B(r)>J  -  2w8(wM°fl°  +  }  _  exp(-/CTwTF{<U^(0),A^(/)1>>  ’  (C4) 

where  0  denotes  commutation.  Using  (C2)  and  (C3)  we  can  write 

lA^<0),A*</)]  -  Z#  |f  [qy(0).q*(/)J  .  (C5) 

J,k"j  "k 

In  die  harmonic  approximation,  qk(t)  in  the  Heisenberg  represenution  is  the  same  as  the  clas¬ 
sical  solution  for  a  simple  harmonic  oscillator 

q*(r)  -  q*(0)coa vkt  +  ^--  sinmkt ,  (C6) 

where  a  is  the  momentum  corresponding  to  the  kth  normal  mode.  Using  this,  we  obtain 
W0),*<r)l  -  8*^-sin«*t  .  (C7) 

Putting  this  back  in  Eq.  (C4),  and  using  the  relations 


sinw*/  - 


•a 

2ar8(«t)  —  J  e~‘“'  dt  , 


we  obtain 


f  (</4(0)2Mf)>)  -  2w(8(«)/l0fl® 


+  1 - t — -  ;  ~  »y)  ~  *(«  +  »>)1 

1  -  exp(-/Jir«v)  j  2 uj  Bqj  dq,  '  1 

For  the  infrared  case,  we  then  have 


(CIO) 


/•|<M.(0)m.(/)>}  -  2it(8(«)(m”)J 


+  1 - 7~  rni  ~  ~  +  » 

1  -  exp(-/37fw)  “  2m j  bqj  >  1 


(Cll) 


where  /*<?  is  the  a  component  of  the  equilibrium  dipole  moment,  and  is  the  angular  fre¬ 
quency  of  the  jth  normal  mode. 

For  the  Raman  case  we  have 


F { <  P„0(O) Pyi(t)  >  |  -  27r(«(a>)/»0a/»°& 


+  *| - 7  ,  Z  T  *1-  —  utj)  —  6(ci*  +  <tfy)]|  , 

1  -  exp(-07M  j  2t»j  dqf  dqf  '  ’  J 

where  P%  is  the  a/3  component  of  the  equilibrium  polarizability  tensor. 


(C12) 
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